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H-Bases  II:  Applications 
to  Numerical  Problems 


H.  Michael  Moller  and  Thomas  Sauer 


Abstract.  We  show  how  H-bases  can  be  applied  to  polynomial  interpo- 
lation and  for  the  solution  of  systems  of  nonlinear  equations.  We  will  give 
an  example  of  a system  of  polynomial  equations  where  the  H-basis  leads 
to  more  stable  computations  than  with  the  Grobner  basis. 


§1.  Introduction 

In  the  preceding  paper  [12],  we  introduced  the  notion  of  H-bases  for  poly- 
nomial ideals,  and  showed  how  to  construct  H-bases  in  the  numerically  most 
interesting  case  of  a zero  dimensional  ideal.  In  this  paper  we  consider  two  prob- 
lems from  Numerical  Analysis,  namely  polynomial  interpolation  and  solving 
systems  of  polynomial  equations,  and  point  out  how  H-bases  can  be  applied 
to  both.  More  precisely,  in  both  cases  the  computation  of  normal  forms  with 
respect  to  an  ideal  plays  a crucial  role,  and  with  the  basic  results  from  [12] 
available,  H-bases  yield  a perfect  replacement  for  the  Grobner  bases  which  are 
normally  and  frequently  used  to  do  this  job  [8] . Finally,  we  will  consider  an 
example  where  a properly  chosen  H-basis  leads  to  a significant  stabilization 
of  the  computations  in  comparison  with  the  use  of  Grobner  bases. 


§2.  Interpolation 

A finite  set  0 C n'  of  linearly  independent  functionals  on  n is  said  to  define 
an  ideal  interpolation  scheme  if  its  kernel,  ker©  C n,  is  an  ideal  in  n.  Given 
an  ideal  interpolation  scheme  © and  a polynomial  f 6 U,  the  interpolation 
problem  consists  of  finding  p € n such  that 

0(p)  = ©(/),  i.e.,  tf(p)  = 1 9(f),  i9  6 0.  (1) 
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So  far,  we  have  put  no  restrictions  on  p;  hence,  there  are  infinitely  many 
solutions  to  (1).  More  precisely,  if  p is  any  solution  of  (1),  hence  f — p £ ker  0, 
then  the  set  of  all  solutions  is  the  equivalence  class 

[p]  = p + ker  0 = / + ker  0 [/] . 

We  denote  the  linear  space  of  all  equivalence  classes  by  If/ ker  0,  and  remark 
that  (dim  11/  ker  0)  = #0.  Of  course,  in  order  to  compute  interpolation  poly- 
nomials, we  must  find  a way  to  choose  a specific  element  from  the  equivalence 
class  [/].  A “natural”  choice  is  to  take  the  normal  form  NF  (/,  H),  where  H 
is  an  H-basis  for  ker0.  Since  [/]  = [g]  implies  that  / — g £ (Tt) , and  since 
NF  (• ,Tt ) is  a linear  operator,  we  have  that 

[/]  = [<?]  =*>  NF  (f,H)  = NF  {g,H)  + W - NF  (g,H) . 

K *>/ * 

=0 

Hence,  NF  ([/] ,Tt)  = NF  (/,  ft),  that  is,  the  normal  form  is  the  same  for  any 
element  of  the  same  equivalence  class.  This  algebraic  approach  also  allows 
for  interpolation  of  functionals  which  are  only  given  implicitly , that  is,  by  an 
ideal  IcH:  compute  an  H-basis  ft  for  X and  the  interpolation  operator  is  the 
“remainder  of  division”  NF  (-,ft).  It  is  worthwhile  to  remark  that  one  of  the 
oldest  papers  on  multivariate  interpolation,  namely  [6],  starts  with  implicitly 
given  interpolation  nodes. 

Another  approach  is  to  look  for  a polynomial  space  Pen  which  allows 
for  unique  interpolation  with  respect  to  0;  to  restrict  the  number  of  solutions 
to  this  problem,  one  usually  demands  the  interpolation  operator  L-p  : n — > V 
to  be  degree  reducing  [3],  that  is, 

deg  Lrf  < deg  /,  / £ n. 

Such  an  interpolation  space  with  a degree  reducing  interpolation  operator 
is  called  a minimal  degree  interpolation  space.  The  most  prominent  minimal 
degree  interpolation  spaces  is  the  least  interpolation  space  introduced  by  de 
Boor  et  al  in  [2],  and  is  the  unique  degree  reducing  interpolation  space  which 
satisfies  the  additional  condition 

V=  fl  ker  <?(£), 

qe ker©  K°Xl  °Xn ' 

On  the  other  hand,  it  is  obvious  that  the  operator  NF  (•,  ft)  is  degree  reducing, 
linear  and  interpolating,  hence  all  the  spaces  V = NF  (n,  ft),  for  any  H- 
basis  ft,  are  minimal  degree  interpolation  spaces  with  interpolation  operator 
Lp  = NF  (•,  H).  Moreover,  it  is  even  possible  to  recover  known  minimal  degree 
interpolation  spaces  by  this  algebraic  process. 

Theorem  1.  [15]  The  least  interpolation  space  is  given  as  NF  where 

Ti  is  an  orthogonal  H-basis  with  respect  to  the  inner-product, 


(p,q)  = (p(D)q)(  o),  p,gen. 
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§3.  Polynomial  System  Solving 

Probably  the  best-known  and  most  frequent  use  of  Grobner  bases  is  for  solv- 
ing polynomial  systems  of  equations,  where  they  form  a core  part  of  literally 
all  available  computer  algebra  systems.  These  systems  of  equations  arise  natu- 
rally in  a geometric  context,  such  as  finding  solutions  of  geometric  constraints 
(for  example,  any  Euclidean  distance  constraint  yields  a quadratic  equation) 
or  “simply”  computing  the  intersection  of  algebraic  curves/surfaces  given  in 
implicit  form.  So,  given  any  finite  set  X € II  one  wants  to  find  the  associated 
algebraic  variety  X 6 IK  (some  algebraic  closure  of  our  underlying  field  K) 
such  that 

X(X)  = 0,  (2) 


that  is, 


fix)  = 0,  xex,  f ex. 


Note  that  the  emphasis  here  is  not  on  finding  one  solution  (which  could,  at 
least  in  the  case  that  ffX  = n,  be  done  by  a Newton  method),  but  on  finding 
all  solutions  and  obtaining  structural  information  about  the  variety.  It  is  easy 
to  see  that  the  variety  is  not  a property  of  the  specific  set  T , but  of  the  ideal 
(X): 

X(X)  = 0 4=>  (X)  (X)  = 0. 


Therefore,  it  may  be  helpful  to  find  particular  bases  for  (X)  which  allow  for 
an  efficient  solution  of  (2).  The  “classical”  implementation  in  most  Computer 
Algebra  systems  relies  on  the  computation  of  elimination  ideals,  which  means 
the  computation  of  a basis  for  the  subideals 


(X)k  = , k = 1,  ...,n, 

where  {X)n  = (X).  In  fact,  this  corresponds  to  transforming  the  original 
problem  X(X)  = 0 into  a triangular  system 

9i(  zi  ) = 0, 

32(  Xi,  x2  ) = 0, 

: (3) 
9m(  Xli  X2 , ■ Xn  ) — 0. 

Once  such  a triangular  system  is  available,  the  solution  strategy  is  obvious: 
determine  the  zeros  of  the  univariate  polynomial  g\  (aq)  and  substitute  them 
into  </2  (',  X2)  which  is  now,  for  for  any  such  zero,  again  a univariate  polynomial 
in  X2,  and  go  on  with  this  procedure.  Moreover,  such  a triangular  basis  can 
indeed  be  computed:  Qiex,  the  reduced  Grobner  basis  for  (X)  with  respect  to 
the  lexicographical  term  order  where  xi  -<  X2  -<  • ■ ■ -<  x„  has  the  property 
that 

Qk  = Q n K [*1, . . . , a*]  c (X)k 

is  a Grobner  basis  for  (X)k  (cf.  [4,  p.  114]). 
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However,  as  nice  as  this  idea  of  successive  elimination  of  variables  sounds, 
there  are  numerous  drawbacks: 


(i)  The  complexity  of  computing  a lexicographical  Grobner  basis  is  tremen- 
dous, and  even  relatively  “simple”  problems  still  exceed  the  limitations 
of  existing  computing  facilities. 

(ii)  There  are  often  several  polynomials  in  a certain  number  of  variables,  that 
is,  the  system  is  not  as  triangular  as  one  would  want  it  to  be. 

(iii)  The  degree  of  the  polynomial  gi  is  usually  very  high.  This  makes  it 
impossible  to  compute  its  zeros  exactly. 

(iv)  The  tempting  idea  to  find  gi’s  zeros  approximately  and  substitute  these 
values  will  not  lead  very  far  since  it  is  well-known  that  the  zeros  of  a 
polynomial  are  usually  quite  ill-conditioned  with  respect  to  its  coefficients 
(cf.  [5,17]). 


So,  the  summary  is  fairly  disappointing:  elimination  methods  do  not  provide 
a good  tool  to  tackle  polynomial  systems  of  equations.  In  particular,  they  rely 
too  much  on  symbolic  methods  (with  exact  computations)  to  become  a useful 
tool  in  numerical  applications. 

A different  approach  has  been  proposed  quite  recently  by  Stetter  [16]  (see 
also  [10];  in  [7]  this  method  is  partly  attributed  to  Stickelberger)  which  is  based 
on  transforming  the  nonlinear  system  of  equations  into  an  eigenvalue  problem 
for  which  a huge  library  of  powerful  routines  is  available.  For  that  purpose, 
let  us  assume  that  the  set  of  solutions  X is  finite  (that  is,  the  associated  ideal 
(IF}  is  zero  dimensional)  and  that  all  the  common  zeros  are  simple.  The  latter 
restriction  is  made  to  keep  the  presentation  simple;  details  on  how  to  handle 
multiplicities  can  be  found  in  [10].  We  first  note  that  for  any  / £ n,  the 
mapping 


*/: 


fn  /(F) 
I [p] 


n /(T) 

[f -p) 


is  a homomorphism  on  the  #X-dimensional  linear  space  n/  (F).  Now,  sup- 
pose for  a moment  that  we  know  X.  Then  there  are  polynomials  px  £ n, 
x € X,  defined  by 

Px  (x  ) ^ X,Xr  , X,  X £ X J 


which  form  a basis  for  11/  (F),  i.e., 


n/  (F)  = span  { [px]  : x £ X } . 

Obviously,  for  any  x £ X,  the  polynomial  gx  — (/  — }(x))px  satisfies  gx(X)  = 
0,  and  therefore 


[0]  = [Px]  = [ (/  - }{x))Px  } = $/  bx ] - f(x)  [Px]  • 


What  we  have  proved  with  this  simple  argument  is  the  following  crucial  the- 


orem. 
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Theorem  2.  The  polynomials  px,  x € X,  are  joint  eigenvectors  of  all  homo- 
morphisms  $f,  / € II,  with  respect  to  the  eigenvalue  f(x). 

This  result  again  suggests  a strategy  to  solve  polynomial  systems  of  equa- 
tions: compute  a set  of  representers  for  11/  (IF},  that  is,  a finite  set  V C II  of 
linearly  independent  polynomials  such  that 

11/  (T)  = span  { [p]  : p € V } , 

and  compute  the  matrix  Mf  which  describes  the  action  of  with  respect  to 
the  basis  V.  The  eigenvalues  of  such  matrices  yield,  when  combined  appro- 
priately, the  solutions  X.  We  remark  that  the  (transpose  of  the)  matrix  Mf 
is  called  the  multiplication  table  for  / with  respect  to  V,  and  that  the  original 
goal  for  Buchberger’s  doctoral  thesis  (supervised  by  Grobner)  was  not  the 
invention  of  Grobner  bases  but  the  computation  of  multiplication  tables.  Of 
course,  the  most  natural  approach  would  be  to  compute  the  multiplication 
tables  MXj , j = 1 ,n,  for  the  coordinate  functions  and  thus  compute  the 
respective  coordinates  of  the  elements  of  X as  the  eigenvalues  of  the  multipli- 
cation table.  Note  that  the  different  components  are  finally  “glued  together” 
by  the  requirement  that  they  must  correspond  to  the  same  eigenvector. 

What  we  now  have  is  the  possibility  of  reducing  the  search  for  the  solu- 
tions of  a polynomial  system  of  equations  to  an  eigenvalue  problem,  provided 
that  we  are  able  to  perform  two  operations: 

(i)  Given  a basis  T for  an  ideal  (T)  compute  a basis  V of  representers  for 

n/(x). 

(ii)  Having  this  basis  available  and  given  any  / € n,  compute  the  multipli- 
cation table  Mf  with  respect  to  V. 

Fortunately,  this  is  where  [12]  enters  - the  answer  are  normal  forms:  if  H is 
an  H-basis  for  (X),  then  any  basis  for  NF  (n,  IT)  is  exactly  the  desired  V,  and 
the  action  of  can  be  computed  by  expanding  NF  (/  • p,  TL)  for  all  p G V, 
which  yields  the  multiplication  table  M f . The  remaining  question  is  “why  H- 
bases?”,  and  this  question  is  justified  since  the  computation  of  normal  forms 
and  thus  of  multiplication  tables  is  perfectly  possible  with  the  help  of  Grobner 
bases  as  well.  To  give  a partial  answer  to  this  question,  we  look  at  an  example. 

§4.  When  Two  Ellipses  Meet 

In  this  section  we  consider  a simple  example  which  will  show  that  also  the 
eigenvalue  method  can  encounter  serious  obstacles,  in  particular  when  Grobner 
bases  are  involved.  The  important  thing  here  is  simplicity,  as  it  will  not  be 
too  surprising  if  extremely  complicated  and  difficult  examples  cause  problems. 

We  consider  the  two  ellipses 

/ (xt  y)  = y 2 - 1, 

g{x,y)  = jjz2  + ^y2  - l. 


338 


H.  M.  Moller  and  T.  Sauer 


Fig.  1.  The  two  ellipses. 


Clearly,  these  two  ellipses  intersect  in  the  four  well-separated  points  (±1,  ±1) 
as  can  be  seen  in  Fig.  1. 

Now,  we  are  going  to  perturb  g a little  bit  and  replace  it  by  g$  = 
g (Aj,(x,y)),  where  A </,  denotes  the  rotation 


A,p(x,  y) 


COS  (f> 

sintf) 

X 

— sin  tf> 

cos  <f> 

y _ 

Note  that  we  have  in  mind  small  values  of  <p,  so  the  intersections  should  still 
be  close  to  (±1,±1)  and  the  problem  should  still  be  well-conditioned. 

Recalling  that  lexicographic  Grobner  bases  are  known  as  troublemakers, 
we  first  try  some  “better”  Grobner  basis,  namely  the  one  which  is  based  on 
the  graded  lexicographic  term  order  with  x -<  y.  Note  that  this  ideal  basis  is 
not  only  a Grobner  basis,  but  also  an  H-basis.  In  this  case  the  Grobner  bases 
G<j,  consists,  for  0 / 0,  of  the  three  polynomials 

4 sin  (f>  xy  + 3 cos  cj>  x2  — 3 cos  <j>, 

x2  + 2 y2  - 3, 

cos  (f>  (cos2  (j>  + 8)  x3  — 3 cos  <$  (cos2  cf>  + 2)  x + 12  sin  cf>  (sin2  <f>  — l)  y, 


while 

Go  = {x2  - l,y2  - 1}  . 

Here  we  already  observe  that  some  singularity  must  appear  for  (f>  = 0,  since 
Go  is  not  just  a limit  cf>  — > 0 of  G<t>,  although  the  basis  changes  continuously 
with  respect  to  (f>.  The  singularity  becomes  more  apparent  if  we  look  at  the 
normal  forms,  which  are 


•p  _ f {l ,x,y,x2}  if  </>  # 0, 
* 1 {l,x,y,xy}  if  ^6  = 0. 
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Finally,  the  multiplication  tables  Mx ^ for  the  multiplication  by  x take  the 
form 


while  the  multiplication  table 


; 

■o 

1 

0 

0- 

ii 

o 

af 

1 

0 

0 

0 

0 

0 

0 

1 

l 

table 

.0 

0 

1 

0. 

■o 

1 

0 

0 


0 

0 

0 

1 


3 cos  <f> 

4 sin  0 

o 

0 

_ 3 cos  <f> 
4 sin  0 


o 

0 5 cos2  0— 8 
cos2  0+8 

1 o sin  tft  cos 

cos2  $+8 

0 


<t>  ¥=  °> 


provides  us  with  difficulties.  Not  only  does  this  matrix  not  converge  to  MXto 
for  <j>  —>  0,  but  some  entries  in  this  matrix  even  diverge  to  ±oo,  respectively. 
Indeed,  if  one  tries  to  compute  the  eigenvalues  and  eigenvectors  of  this  matrix 
for  small  values  of  0,  things  become  disastrous:  A Maple  computation  with 
10  digits  worked  until  about  0 ~ 10— 5 , where  an  error  message  reported  that 
the  QR  algorithm  did  not  work.  For  smaller  values,  like  (j>  ~ 10~6,  Maple 
invented  complex  zeros  with  an  imaginary  part  of  the  magnitude  0.5  x 10~5 
which  by  far  exceeds  any  negligible  machine  number.  On  the  other  hand, 
Octave,  a free  Matlab  clone  whose  Linear  Algebra  facilities  are  based  on 
LAPACK  [1],  reproduced  the  eigenvalues  correctly,  but  gave  eigenvectors  which 
were  practically  0. 

Hence,  we  end  up  with  some  kind  of  paradox  which  is  due  to  a singularity 
at  cf)  = 0:  though  the  original  problem  of  solving  the  polynomial  system  of 
equations  is  very  well-conditioned,  the  graded  lexicographical  Grobner  basis 
is  extremely  sensitive  to  very  small  perturbations  (|0|  < HR5),  but  by  far  not 
so  sensitive  to  relatively  “large”  (|0[  > 10-5)  perturbations. 

Similar  problems  appear  when  we  replace  the  graded  lexicographical 
Grobner  basis  by  a purely  lexicographical  one  with  x -<  y which  yields  the 
normal  forms 

■p  _ f {l,x,x2,x3}  if  0^0, 

0 \ {l,x,y,xy}  if  0 = 0. 

Though  the  components  of  the  multiplication  table  Mx^  at  least  are  contin- 
uous functions  in  0 and  remain  bounded  in  this  case,  the  limit  0 — > 0 again  is 
not  MXfi.  But  the  multiplication  tables  My^  with  respect  to  the  purely  lexi- 
cographical Grobner  basis  is  even  worse:  its  entries  are  either  zero  or  diverge 
for  0 — > 0. 

The  behavior  of  the  Grobner  bases  at  0 = 0 raises  the  question  of  whether 
this  singularity  is  systematic,  that  is,  intrinsic  to  the  problem,  or  if  it  is  a repre- 
sentation singularity  generated  by  the  Grobner  bases.  Systematic  singularities 
appear,  for  example,  if  several  zeros  “collapse”  into  one  multiple  zero  which 
leads  to  extremely  intricate  problems  in  the  multivariate  case  [9].  Here,  how- 
ever, the  good  separation  of  the  zeros  suggests  the  conjecture  that  we  only 
face  a representation  singularity. 
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Indeed,  since  H-bases  leave  more  degrees  of  freedom,  we  can  try  another 
one  which  is  now  based  on  orthogonalization.  For  this  purpose,  we  use  the 
inner-product 

(P.9)  = (p{D)  q ) (0) 

and  recall  from  [11,  Theorem  5.3]  that  the  set  {/,  g#}  is  already  an  H-basis. 
Moreover,  the  normal  form  space,  which  is,  according  to  Theorem  1,  the  least 
interpolation  space,  is  spanned  by 

= { 1,  x,  y,  2 sirup  x 2 — 3 cos  <f>  xy  — sin  <f>  y2} 


and  depends  continuously  on  $ with 


lim  7^  = Vo  = {1  ,x,y,xy}. 

<p — >0 


Then  one  can  compute  the  respective  multiplication  table  as 


•o 

1 + £l((p) 

£2(0) 

£3(0)  ' 

1 

0 

0 

£4(0) 

1 

0 

0 

1 + £5(0) 

.0 

eeM 

1 + £7  W 

E8(lp)  . 

where  Sj  (■),  j = 1, . . . , 8,  are  continuous  functions  which  vanish  at  the  origin. 
In  particular,  M*  ^ — > MXi o as  a:  — ► 0 and  the  computation  of  eigenvalues 
and  eigenvectors  of  ^ can  now  be  done  with  sufficient  accuracy.  However, 
we  remark  that  the  fact  that  the  matrices  M*  ^ and  M*  ^ have  two  approxi- 
mately double  eigenvalues  ±1,  requires  some  extra  care  when  connecting  these 
individual  values  in  the  final  determination  of  the  intersections. 


§5.  Summary 

We  have  given  examples  of  numerical  applications  which  can  be  reduced  to 
the  computation  of  normal  forms  with  respect  to  a certain  polynomial  ideal, 
an  operation  which  is  usually  performed  using  a Grobner  basis.  On  the  other 
hand,  H-bases  could  as  well  be  used  for  normal  form  computations,  and  their 
greater  flexibility  may  yield  stabilizing  effects  which  are  highly  desired  in  nu- 
merical computations. 
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